# Project Information---
#  Project: Educating for Inclusion: Diversity education programs can reduce prejudice towards outgroups in Israel
#  Journal: PNAS
#  Last updated: Tue Mar  7 15:42:43 2023
#  Purpose: Estimate all results for paper -- Study I
#  Outputs: All figures and tables associated with Study 1 in paper and appendix
#  Machine: Chagai Weiss's Macbook Pro


#Load Packages------
library("tidyverse")
library("pmdplyr")
library("DataCombine")
library("dplyr")
library("readxl")
library("lubridate")
library("estimatr")
library("ltm")
library("effectsize")
library("ggthemes")
library("stargazer")
library("fwildclusterboot")
library("RItools")
library("xtable")
library("multiwayvcov")
library("lmtest")
library("RColorBrewer")
library("interflex")
library("ri2")
library("exr")

set.seed(999)

# Read Data------
setwd(paste(getwd()))
full_data <- read.csv("study1_data.csv")

# Results in Main Text-------
## Figure 3: Main ATE Study 1------

###WCB Thermometer-----
therm_lm <- 
 full_data %>% 
 mutate(.,
        therm_index_post = standardize(therm_index_post)) %>% 
 lm(therm_index_post ~ treatment + therm_index_pre +
     na_therm_arab + na_therm_imgrnt + na_therm_uo + na_therm_blind+
     grade.x + boy.x,
    data = .) 

therm_bstrp <- 
 cluster.boot(therm_lm, full_data$class.x, 
              parallel = T, R = 1000, 
              boot_type = "wild")

therm_bstrp_out <- 
 coeftest(therm_lm, therm_bstrp) %>% 
 tidy(.,
      conf.int = T) %>% 
 mutate(.,
        Outcome = "Thermometers") %>% 
 filter(.,
        term == "treatment")

###WCB Similarity-----
sim_lm <- 
 full_data %>% 
 mutate(.,
        similarity_index_post = standardize(similarity_index_post)) %>% 
 lm(similarity_index_post ~ treatment + similarity_index_pre + 
     na_similar_blind + na_similar_imgrnt + na_similar_arab + na_similar_uo+
     grade.x + boy.x,
    data = .) 

sim_bstrp <- 
 cluster.boot(sim_lm, full_data$class.x, 
              parallel = T, R = 1000, 
              boot_type = "wild")

sim_bstrp_out <- 
 coeftest(sim_lm , sim_bstrp) %>% 
 tidy(.,
      conf.int = T) %>% 
 mutate(.,
        Outcome = "Similarity  \nPerceptions") %>% 
 filter(.,
        term == "treatment")

###WCB Contact Index-----
cntct_ind_lm <- 
 full_data %>% 
 mutate(.,
        contact_general_post = standardize(contact_general_post)) %>% 
 lm(contact_general_post ~ treatment + contact_general_pre + 
     na_cntnct_blin_play + na_cntnct_blind_bday + na_cntnct_blind_hwork + na_cntnct_blind_lost+
     na_cntnct_arab_bday + na_cntnct_arab_hwork + na_cntnct_arab_lost + na_cntnct_arab_play +
     na_cntnct_imgrnt_bday + na_cntnct_imgrnt_hwork + na_cntnct_imgrnt_play + na_cntnct_imgrnt_lost+
     na_cntnct_uo_bday + na_cntnct_uo_hwork + na_cntnct_uo_play + na_cntnct_uo_lost +
    grade.x + boy.x,
    data = .) 

cntct_ind_bstrp <- 
 cluster.boot(cntct_ind_lm, full_data$class.x, 
              parallel = T,  R = 1000, 
              boot_type = "wild")

cnt_ind_bstrp_out <- 
 coeftest(cntct_ind_lm , cntct_ind_bstrp) %>% 
 tidy(.,
      conf.int = T) %>% 
 mutate(.,
        Outcome = "Contact \nIntentions") %>% 
 filter(.,
        term == "treatment")


###WCB Contact Behavior-----

cntct_bh_lm <- 
 full_data %>% 
 mutate(.,
        attend_contact = standardize(attend_contact)) %>% 
 lm(attend_contact ~ treatment + 
     contact_general_pre + therm_index_pre + similarity_index_pre +div_general_scale_pre+
     na_similar_blind + na_similar_imgrnt + na_similar_arab + na_similar_uo+
     na_div_benefit + na_div_enjoy_learn + na_div_learn + na_div_listen + na_div_enjoy_play+
     na_therm_arab + na_therm_blind + na_therm_imgrnt + na_therm_imgrnt+
     na_cntnct_blin_play + na_cntnct_blind_bday + na_cntnct_blind_hwork + na_cntnct_blind_lost+
     na_cntnct_arab_bday + na_cntnct_arab_hwork + na_cntnct_arab_lost + na_cntnct_arab_play +
     na_cntnct_imgrnt_bday + na_cntnct_imgrnt_hwork + na_cntnct_imgrnt_play + na_cntnct_imgrnt_lost+
     na_cntnct_uo_bday + na_cntnct_uo_hwork + na_cntnct_uo_play + na_cntnct_uo_lost +
     grade.x + boy.x,
    data = .) 

cntct_bh_bstrp <- 
 cluster.boot(cntct_bh_lm, full_data$class.x, 
              parallel = T, R = 1000, 
              boot_type = "wild")

cnt_bh_bstrp_out <- 
 coeftest(cntct_bh_lm , cntct_bh_bstrp) %>% 
 tidy(.,
      conf.int = T) %>% 
 mutate(.,
        Outcome = "Register \nContact") %>% 
 filter(.,
        term == "treatment")




###WCB Diversity index-----
div_lm <- 
 full_data %>% 
 mutate(.,
        div_general_scale_post = standardize(div_general_scale_post)) %>% 
 lm(div_general_scale_post ~ treatment+ div_general_scale_pre + 
     na_div_benefit + na_div_enjoy_learn + na_div_learn + na_div_listen + na_div_enjoy_play+ 
   grade.x + boy.x,
    data = .) 

div_bstrp <- 
 cluster.boot(div_lm, full_data$class.x, 
              parallel = T, use_white = NULL, R = 1000, 
              boot_type = "wild")

div_bstrp_out <- 
 coeftest(div_lm , div_bstrp) %>% 
 tidy(.,
      conf.int = T) %>% 
 mutate(.,
        Outcome = "Diversity \nAttitudes") %>% 
 filter(.,
        term == "treatment")

#### Plot Bootsrap Models-----

# Combine all models together
bootstrap_models <- 
 rbind(therm_bstrp_out, sim_bstrp_out, 
       cnt_ind_bstrp_out, cnt_bh_bstrp_out,div_bstrp_out) %>% 
 mutate(.,
        se = paste0("(", round(std.error,2),")"),
        pese = paste(round(estimate,2), se))
# set order
bootstrap_models$Outcome <- 
 factor(bootstrap_models$Outcome, 
        c("Register \nContact", "Similarity  \nPerceptions",
         "Diversity \nAttitudes",
          "Contact \nIntentions",
           "Thermometers"))



ggplot(bootstrap_models, aes(x=estimate, y=Outcome,
                             label = estimate, color = "dodgerblue4")) + 
 geom_errorbar(aes(xmin=conf.low, xmax=conf.high), width=.1) +
 geom_vline(xintercept = 0,
            color = "black",
            linetype = "dashed",
            size = 0.4) +
 geom_point()+
 geom_text(aes(label= pese), 
           position = position_dodge(width = .4),
           hjust = -3, 
           family = "Times",
           size = 3) +
 labs(x = "",
      y = "")+
 xlim(-.2,1)+
 scale_color_manual(values = c("dodgerblue4"))+
 theme(text = element_text(size = 16, family = "Times"),
       legend.position = "none",
       panel.grid.major = element_blank(),
       panel.background = element_blank())

#Appendix------

##Figure A1, Panels a/b: Levels of pre-treatment thermometor/similarity-------
###Therm 
# select relevant measures
therm_data <- full_data %>% 
 dplyr::select(.,
               blind_therm_pre, arab_therm_pre,
               uo_therm_pre, immigrant_therm_pre)

# pivot data
therm_piv <-
 therm_data %>% 
 pivot_longer(blind_therm_pre:immigrant_therm_pre, 
              names_to = "group", 
              values_to = "value") 
# get mean & SD
therm_estimates <- 
 therm_piv %>%
 group_by(group) %>%
 summarize(Mean_therm = mean(value, na.rm=TRUE),
           Sd_therm = sd(value, na.rm=TRUE),
           n_therm = n()) %>% 
 mutate(se_therm = Sd_therm / sqrt(n_therm),
        lower_ci_rust = Mean_therm - qt(1 - (0.05 / 2), n_therm- 1) * se_therm,
        upper_ci_rust = Mean_therm + qt(1 - (0.05 / 2), n_therm - 1) * se_therm) %>% 
 mutate(.,
        group = case_when(
         group == "blind_therm_pre" ~ "Visually \nImpaired",
         group == "arab_therm_pre" ~ "Arab",
         group == "uo_therm_pre" ~ "Ultra-Orthodox",
         group == "immigrant_therm_pre" ~ "Immigrant"
        ))


ggplot(therm_estimates, aes(x=Mean_therm, y=group, colour="dodgerblue4")) + 
 geom_errorbar(aes(xmin=lower_ci_rust, xmax=upper_ci_rust), width=.1) +
 geom_point()+
 labs(x = "Thermometer",
      y = "")+
 geom_vline(xintercept = mean(therm_piv$value, na.rm = T),
            color = "black",
            linetype = "dashed",
            size = 0.4)+
 xlim(0,100)+
 scale_color_manual(values = c("dodgerblue4"))+
 theme(text = element_text(size = 20, family = "Times"),
       legend.position = "none",
       panel.grid.major = element_blank(),
       panel.background = element_blank())


###Similarity
# select relevant measures
sim_data <- full_data %>% 
 dplyr::select(.,
               similarity_arabs_pre, similarity_uo_pre,
               similarity_blind_pre, similarity_imgrnt_pre)

# pivot data
sim_piv <-
 sim_data %>% 
 pivot_longer(similarity_arabs_pre:similarity_imgrnt_pre, 
              names_to = "group", 
              values_to = "value") 
# get mean & SD
sim_estimates <- 
 sim_piv %>%
 group_by(group) %>%
 summarize(Mean_sim = mean(value, na.rm=TRUE),
           Sd_sim = sd(value, na.rm=TRUE),
           n_sim = n()) %>% 
 mutate(se_sim = Sd_sim / sqrt(n_sim),
        lower_ci_rust = Mean_sim - qt(1 - (0.05 / 2), n_sim- 1) * se_sim,
        upper_ci_rust = Mean_sim + qt(1 - (0.05 / 2), n_sim - 1) * se_sim) %>% 
 mutate(.,
        group = case_when(
         group == "similarity_blind_pre" ~ "Visually \nImpaired",
         group == "similarity_arabs_pre" ~ "Arab",
         group == "similarity_uo_pre" ~ "Ultra-Orthodox",
         group == "similarity_imgrnt_pre" ~ "Immigrant"
        ))


ggplot(sim_estimates, aes(x=Mean_sim, y=group, colour="dodgerblue4")) + 
 geom_errorbar(aes(xmin=lower_ci_rust, xmax=upper_ci_rust), width=.1) +
 geom_point()+
 labs(x = "Similarity Index",
      y = "")+
 scale_color_manual(values = c("dodgerblue4"))+
 geom_vline(xintercept = mean(sim_piv$value, na.rm = T),
            color = "black",
            linetype = "dashed",
            size = 0.4)+
 xlim(1,5)+
 theme(text = element_text(size = 20, family = "Times"),
       legend.position = "none",
       panel.grid.major = element_blank(),
       panel.background = element_blank())


##Figure A4: Plot of days between treatment and outcome collection-----

full_data$days_since_intervention <- as.numeric(difftime(as.Date(full_data$StartDate.y),
                                                         as.Date("2021-06-03"),
                                                         units=c("days")))
mean <- mean(full_data$days_since_intervention[full_data$treatment==1], na.rm = T)
full_data %>% 
 filter(.,
        treatment == 1) %>% 
 ggplot(., aes(x=days_since_intervention, fill = "dodgerblue4")) +
 geom_histogram(alpha = .6)+
 geom_vline(aes(xintercept=mean), color = "dodgerblue4",
            linetype="dashed")+
 scale_x_continuous(breaks = scales::pretty_breaks(n = 10))+
 labs(x= "Days Since Treatment",
      y = "Count",
      color = "",
      fill = "")+
 scale_fill_manual(values = c("dodgerblue4"))+
 scale_color_manual(values = c("dodgerblue3"))+
 theme(panel.spacing = unit(1.5, "lines"),
       text = element_text(size = 14, family = "Times"),
       # axis.text.y = element_text(size=14),
       legend.position = "none",
       panel.grid.major = element_blank(),
       panel.background = element_blank(),
       legend.key=element_blank())


##Table A2: Descriptive Statistics-------
# create descriptive stats dataframe
disc_table <- full_data %>% 
 mutate(.,
        grade_4 = ifelse(grade.x == 4,1,0),
        grade_5 = ifelse(grade.x == 5,1,0),
        grade_6 = ifelse(grade.x == 6,1,0)) %>% 
 dplyr::select(.,
               boy.x, age.x, grade_4,grade_5, grade_6)

# Output table
stargazer(as.data.frame(disc_table),
          covariate.labels = c("Boy", "Age", "Grade 4",
                               "Grade 5", "Grade 6"),
          title = "Descriptive Statistics",
          label = "tab:dsc_experiment1",
          style = "qje",
          #  notes = c("Therm refers to feelings thermometer, SD refers to social disctance scale, Manip refers to manipulation check."),
          notes.append = T,
          notes.align = "l",
          omit.summary.stat = c("p25", "p75"))



##Table A3: Balance Check------
 balanceTest(treatment~ boy.x + age.x+
              therm_index_pre +
              arab_therm_pre +immigrant_therm_pre + uo_therm_pre + blind_therm_pre +
              similarity_index_pre+
              similarity_uo_pre + similarity_arabs_pre + similarity_imgrnt_pre + similarity_blind_pre + 
              div_general_scale_pre +contact_general_pre + cluster(class.x),
             data = full_data) 

##Figure A11: Correlates of non-response------

# Create missingness indicators 
full_data <-
 full_data %>% 
 mutate(.,
        no_uo_contact = ifelse(cntct_uo_answers ==0,1,0),
        no_arab_contact = ifelse(cntct_arab_answers ==0,1,0),
        no_blind_contact = ifelse(cntct_blind_answers ==0,1,0),
        no_imgr_contact = ifelse(cntct_imgrnt_answers ==0,1,0),
        no_cntct_indx = ifelse(is.na(contact_general_post),1,0),
        no_therm_post = ifelse(is.na(therm_index_post),1,0),
        no_sim_post = ifelse(is.na(similarity_index_post),1,0),
        no_div_post = ifelse(is.na(div_general_scale_post),1,0)
 )

# Model atrition
no_cntct_blind <- 
 lm(no_blind_contact ~
     treatment + boy.x + age.x + contact_general_pre, 
    data = full_data)
no_cntct_blind_out <-
 tidy(no_cntct_blind, conf.int = T) %>% 
 filter(.,
        term %in% c("treatment", "contact_general_pre")) %>% 
 mutate(.,
        Outcome = "Blind Contact",
        term = case_when(
         term == "treatment" ~ "Treatment",
         term == "contact_general_pre" ~ "Pre-Treatment \nOutcome",
        ))

no_cntct_arab <- 
 lm(no_arab_contact ~
     treatment + boy.x + age.x + contact_general_pre, 
    data = full_data)
no_cntct_arab_out <-
 tidy(no_cntct_arab, conf.int = T) %>% 
 filter(.,
        term %in% c("treatment", "contact_general_pre")) %>% 
 mutate(.,
        Outcome = "Arab Contact",
        term = case_when(
         term == "treatment" ~ "Treatment",
         term == "contact_general_pre" ~ "Pre-Treatment \nOutcome",
        ))

no_cntct_imgrnt <- 
 lm(no_imgr_contact ~
     treatment + boy.x + age.x + contact_general_pre, 
    data = full_data)
no_cntct_imgrnt_out <-
 tidy(no_cntct_imgrnt, conf.int = T) %>% 
 filter(.,
        term %in% c("treatment", "contact_general_pre")) %>% 
 mutate(.,
        Outcome = "Immigrant Contact",
        term = case_when(
         term == "treatment" ~ "Treatment",
         term == "contact_general_pre" ~ "Pre-Treatment \nOutcome",
        ))

no_cntct_uo <- 
 lm(no_uo_contact ~
     treatment + boy.x + age.x + contact_general_pre, 
    data = full_data)
no_cntct_uo_out <-
 tidy(no_cntct_uo, conf.int = T) %>% 
 filter(.,
        term %in% c("treatment", "contact_general_pre")) %>% 
 mutate(.,
        Outcome = "UO Contact",
        term = case_when(
         term == "treatment" ~ "Treatment",
         term == "contact_general_pre" ~ "Pre-Treatment \nOutcome",
        ))


no_cntct_indx <- 
 lm(no_cntct_indx ~
     treatment + boy.x + age.x + contact_general_pre, 
    data = full_data)
no_cntct_indx_out <-
 tidy(no_cntct_indx, conf.int = T) %>% 
 filter(.,
        term %in% c("treatment", "contact_general_pre")) %>% 
 mutate(.,
        Outcome = "Contact Index",
        term = case_when(
         term == "treatment" ~ "Treatment",
         term == "contact_general_pre" ~ "Pre-Treatment \nOutcome",
        ))


no_therm <-
 lm(no_therm_post ~
     treatment + boy.x + age.x + therm_index_pre, 
    data = full_data)
no_therm_out <-
 tidy(no_therm, conf.int = T) %>% 
 filter(.,
        term %in% c("treatment", "therm_index_pre")) %>% 
 mutate(.,
        Outcome = "Thermometer Index",
        term = case_when(
         term == "treatment" ~ "Treatment",
         term == "therm_index_pre" ~ "Pre-Treatment \nOutcome",
        ))

no_sim <-
 lm(no_sim_post ~
     treatment + boy.x + age.x + similarity_index_pre, 
    data = full_data)
no_sim_out <-
 tidy(no_sim, conf.int = T) %>% 
 filter(.,
        term %in% c("treatment", "similarity_index_pre")) %>% 
 mutate(.,
        Outcome = "Similarity Index",
        term = case_when(
         term == "treatment" ~ "Treatment",
         term == "similarity_index_pre" ~ "Pre-Treatment \nOutcome",
        ))

no_div <-
 lm(no_div_post ~
     treatment + boy.x + age.x + div_general_scale_pre, 
    data = full_data)
no_div_out <-
 tidy(no_div, conf.int = T) %>% 
 filter(.,
        term %in% c("treatment", "div_general_scale_pre")) %>% 
 mutate(.,
        Outcome = "Diversity Index",
        term = case_when(
         term == "treatment" ~ "Treatment",
         term == "div_general_scale_pre" ~ "Pre-Treatment \nOutcome",
        ))

att_models <- 
 rbind(no_cntct_blind_out, no_cntct_arab_out, no_cntct_imgrnt_out,
       no_cntct_uo_out,no_therm_out, no_sim_out, no_div_out, no_cntct_indx_out) %>%
 mutate(.,
        term = case_when(
         term == "Treatment" ~ "Correlation of \nTreatment with \nNonresponse",
         term == "Pre-Treatment \nOutcome" ~ "Correlation of \nPre-Treatment \nOutcome with Nonresponse"
        ))
att_models$term <- 
 factor(att_models$term,
        c("Correlation of \nTreatment with \nNonresponse", "Correlation of \nPre-Treatment \nOutcome with Nonresponse"))
att_order<-
 c("Thermometer Index", 
   "Similarity Index", "Diversity Index",
   "Contact Index",
   "Arab Contact", "Blind Contact",
   "UO Contact", "Immigrant Contact") %>% rev()
ggplot(att_models, aes(x=estimate, y=Outcome, color = "dodgerblue4")) + 
 geom_point(position = position_dodge(width = 0.6), show.legend = FALSE)+
 geom_errorbar(aes(xmin=conf.low, xmax=conf.high), 
               width=.1,
               position = position_dodge(width = 0.6),
               show.legend = FALSE) +
 geom_vline(xintercept = 0,
            color = "black",
            linetype = "dashed",
            size = 0.4) +
 xlim(-.3,.3)+
 scale_y_discrete(limits = att_order)+
 labs(x = "",
      y = "",
      color = "")+
 scale_color_manual(values = c("dodgerblue4"))+
 facet_grid(cols = vars(term))+
 # theme_tufte()+
 #guides(color = guide_legend(reverse = TRUE))+
 theme(panel.spacing = unit(1.5, "lines"),
       text = element_text(size = 12, family = "Times"),
       # axis.text.y = element_text(size=14),
       legend.position = "bottom",
       panel.grid.major = element_blank(),
       panel.background = element_blank(),
       legend.key=element_blank())

##Table A4: Randomization Inference-----
###Estimate RI-----
#Thermometer
# Omit Outcome NAs and rename treatment and outcome
therm_ri_data <-
 full_data %>% 
 filter(.,
        !is.na(therm_index_post)) %>% 
 mutate(.,
        Y = standardize(therm_index_post),
        Z = treatment)

# Declare assignment
declaration_therm <- 
 with(therm_ri_data,{
  declare_ra(
   blocks = grade.x,
   clusters = class.x)
 })

# Run RI
ri_therm <- 
 conduct_ri(
  Y ~ Z + therm_index_pre + boy.x + grade.x +
   na_therm_arab + na_therm_imgrnt + na_therm_uo + na_therm_blind,
  sharp_hypothesis = 0,
  declaration = declaration_therm,
  data = therm_ri_data)
ri_therm_out <- tidy(ri_therm) %>% 
 mutate(.,
        term = "Thermometers")




#Similarity
# Omit Outcome NAs and rename treatment and outcome
sim_ri_data <-
 full_data %>% 
 filter(.,
        !is.na(similarity_index_post)) %>% 
 mutate(.,
        Y = standardize(similarity_index_post),
        Z = treatment)

# Declare assingment
declaration_sim <- 
 with(sim_ri_data,{
  declare_ra(
   blocks = grade.x,
   clusters = class.x)
 })

# Run RI
ri_sim <- 
 conduct_ri(
  Y ~ Z + similarity_index_pre + boy.x + grade.x +
   na_similar_blind + na_similar_imgrnt + na_similar_arab + na_similar_uo,
  sharp_hypothesis = 0,
  declaration = declaration_sim,
  data = sim_ri_data)
ri_sim_out <- tidy(ri_sim) %>% 
 mutate(.,
        term = "Similarity Perceptions")



#Contact Index
# Omit Outcome NAs and rename treatment and outcome
cntct_indx_ri_data <-
 full_data %>% 
 filter(.,
        !is.na(contact_general_post)) %>% 
 mutate(.,
        Y = standardize(contact_general_post),
        Z = treatment)

# Declare assignment
declaration_cntct <- 
 with(cntct_indx_ri_data,{
  declare_ra(
   blocks = grade.x,
   clusters = class.x)
 })

# Run RI
ri_cntct <- 
 conduct_ri(
  Y ~ Z + contact_general_pre + boy.x + grade.x +
  na_cntnct_blin_play + na_cntnct_blind_bday + na_cntnct_blind_hwork + na_cntnct_blind_lost+
   na_cntnct_arab_bday + na_cntnct_arab_hwork + na_cntnct_arab_lost + na_cntnct_arab_play +
   na_cntnct_imgrnt_bday + na_cntnct_imgrnt_hwork + na_cntnct_imgrnt_play + na_cntnct_imgrnt_lost+
   na_cntnct_uo_bday + na_cntnct_uo_hwork + na_cntnct_uo_play + na_cntnct_uo_lost,
  sharp_hypothesis = 0,
  declaration = declaration_cntct,
  data = cntct_indx_ri_data)
ri_cntct_out <- tidy(ri_cntct) %>% 
 mutate(.,
        term = "Contact Intensions")



#Contact Behavior
# Omit Outcome NAs and rename treatment and outcome
cntct_bhv_ri_data <-
 full_data %>% 
 filter(.,
        !is.na(attend_contact)) %>% 
 mutate(.,
        Y = standardize(attend_contact),
        Z = treatment)

# Declare assingment
declaration_cntct_bhv <- 
 with(cntct_bhv_ri_data,{
  declare_ra(
   blocks = grade.x,
   clusters = class.x)
 })

# Run RI
ri_cntct_bhv <- 
 conduct_ri(
  Y ~ Z + contact_general_pre + boy.x + grade.x + 
   therm_index_pre + div_general_scale_pre + similarity_index_pre  +
   na_similar_blind + na_similar_imgrnt + na_similar_arab + na_similar_uo+
   na_div_benefit + na_div_enjoy_learn + na_div_learn + na_div_listen + na_div_enjoy_play+
   na_therm_arab + na_therm_blind + na_therm_imgrnt + na_therm_imgrnt+
   na_cntnct_blin_play + na_cntnct_blind_bday + na_cntnct_blind_hwork + na_cntnct_blind_lost+
   na_cntnct_arab_bday + na_cntnct_arab_hwork + na_cntnct_arab_lost + na_cntnct_arab_play +
   na_cntnct_imgrnt_bday + na_cntnct_imgrnt_hwork + na_cntnct_imgrnt_play + na_cntnct_imgrnt_lost+
   na_cntnct_uo_bday + na_cntnct_uo_hwork + na_cntnct_uo_play + na_cntnct_uo_lost,
  sharp_hypothesis = 0,
  declaration = declaration_cntct_bhv,
  data = cntct_bhv_ri_data)
ri_cntct_bhv_out <- tidy(ri_cntct_bhv) %>% 
 mutate(.,
        term = "Register Contact")


#Diversity Index
# Omit Outcome NAs and rename treatment and outcome
div_ri_data <-
 full_data %>% 
 filter(.,
        !is.na(div_general_scale_post)) %>% 
 mutate(.,
        Y = standardize(div_general_scale_post),
        Z = treatment)

# Declare assingment
declaration_div <- 
 with(div_ri_data,{
  declare_ra(
   blocks = grade.x,
   clusters = class.x)
 })

# Run RI
ri_div <- 
 conduct_ri(
  Y ~ Z + div_general_scale_pre + boy.x + grade.x +
   na_div_benefit + na_div_enjoy_learn + na_div_learn + na_div_listen + na_div_enjoy_play,
  sharp_hypothesis = 0,
  declaration = declaration_div,
  data = div_ri_data)
ri_div_out <- tidy(ri_div) %>% 
 mutate(.,
        term = "Diversity Attitudes")


#### Create Table------

 rbind(ri_therm_out, ri_cntct_out,ri_div_out, ri_sim_out, ri_cntct_bhv_out) %>% 
 rename(.,
        Term = term,
        Estimate = estimate) %>% 
 dplyr::select(.,
               Term:p.value)




##Figure A12: ATE on Disagregate Scales----
###Thermometer-----
####Arabs
therm_arab_lm <- 
 full_data %>% 
 mutate(.,
        arab_therm_post = standardize(arab_therm_post)) %>% 
 lm(arab_therm_post ~ treatment + arab_therm_pre +
     na_therm_arab + grade.x + boy.x,
    data = .) 

therm_arab_bstrp <- 
 cluster.boot(therm_arab_lm, full_data$class.x, 
              parallel = T, R = 1000, 
              boot_type = "wild")

therm_arab_bstrp_out <- 
 coeftest(therm_arab_lm, therm_arab_bstrp) %>% 
 tidy(.,
      conf.int = T) %>% 
 mutate(.,
        Outcome = "Thermometers",
        Group = "Arabs") %>% 
 filter(.,
        term == "treatment")



####Immigrants
therm_imgrnt_lm <- 
 full_data %>% 
 mutate(.,
        immigrant_therm_post = standardize(immigrant_therm_post)) %>% 
 lm(immigrant_therm_post ~ treatment + immigrant_therm_pre +
     na_therm_imgrnt + grade.x + boy.x,
    data = .) 

therm_imgrnt_bstrp <- 
 cluster.boot(therm_imgrnt_lm, full_data$class.x, 
              parallel = T, R = 1000, 
              boot_type = "wild")

therm_imgrnt_bstrp_out <- 
 coeftest(therm_imgrnt_lm, therm_imgrnt_bstrp) %>% 
 tidy(.,
      conf.int = T) %>% 
 mutate(.,
        Outcome = "Thermometers",
        Group = "Immigrants") %>% 
 filter(.,
        term == "treatment")




####UO
therm_uo_lm <- 
 full_data %>% 
 mutate(.,
        uo_therm_post = standardize(uo_therm_post)) %>% 
 lm(uo_therm_post ~ treatment + uo_therm_pre +
     na_therm_uo + grade.x + boy.x,
    data = .) 

therm_uo_bstrp <- 
 cluster.boot(therm_uo_lm, full_data$class.x, 
              parallel = T, R = 1000, 
              boot_type = "wild")

therm_uo_bstrp_out <- 
 coeftest(therm_uo_lm, therm_uo_bstrp) %>% 
 tidy(.,
      conf.int = T) %>% 
 mutate(.,
        Outcome = "Thermometers",
        Group = "Ultra-Orthodox") %>% 
 filter(.,
        term == "treatment")

####Blind
therm_blind_lm <- 
 full_data %>% 
 mutate(.,
        blind_therm_post = standardize(blind_therm_post)) %>% 
 lm(blind_therm_post ~ treatment + blind_therm_pre +
     na_therm_blind + grade.x + boy.x,
    data = .) 

therm_blind_bstrp <- 
 cluster.boot(therm_blind_lm, full_data$class.x, 
              parallel = T, R = 1000, 
              boot_type = "wild")

therm_blind_bstrp_out <- 
 coeftest(therm_blind_lm, therm_blind_bstrp) %>% 
 tidy(.,
      conf.int = T) %>% 
 mutate(.,
        Outcome = "Thermometers",
        Group = "Visually Impaired") %>% 
 filter(.,
        term == "treatment")


###Similarity-----
####Arabs
sim_arabs_lm <- 
 full_data %>% 
 mutate(.,
        similarity_arabs_post = standardize(similarity_arabs_post)) %>% 
 lm(similarity_arabs_post ~ treatment + similarity_arabs_pre + 
     na_similar_arab +
     grade.x + boy.x,
    data = .) 

sim_arab_bstrp <- 
 cluster.boot(sim_arabs_lm, full_data$class.x, 
              parallel = T, R = 1000, 
              boot_type = "wild")

sim_arab_bstrp_out <- 
 coeftest(sim_arabs_lm , sim_arab_bstrp) %>% 
 tidy(.,
      conf.int = T) %>% 
 mutate(.,
        Outcome = "Similarity  \nPerceptions",
        Group = "Arabs") %>% 
 filter(.,
        term == "treatment")



####Immigrants
sim_imgrnt_lm <- 
 full_data %>% 
 mutate(.,
        similarity_imgrnt_post = standardize(similarity_imgrnt_post)) %>% 
 lm(similarity_imgrnt_post ~ treatment + similarity_imgrnt_pre + 
     na_similar_imgrnt +
     grade.x + boy.x,
    data = .) 

sim_imgrnt_bstrp <- 
 cluster.boot(sim_imgrnt_lm, full_data$class.x, 
              parallel = T, R = 1000, 
              boot_type = "wild")

sim_imgrnt_bstrp_out <- 
 coeftest(sim_imgrnt_lm , sim_imgrnt_bstrp) %>% 
 tidy(.,
      conf.int = T) %>% 
 mutate(.,
        Outcome = "Similarity  \nPerceptions",
        Group = "Immigrants") %>% 
 filter(.,
        term == "treatment")



####UO
sim_uo_lm <- 
 full_data %>% 
 mutate(.,
        similarity_uo_post = standardize(similarity_uo_post)) %>% 
 lm(similarity_uo_post ~ treatment + similarity_uo_pre + 
     na_similar_uo +
     grade.x + boy.x,
    data = .) 

sim_uo_bstrp <- 
 cluster.boot(sim_uo_lm, full_data$class.x, 
              parallel = T, R = 1000, 
              boot_type = "wild")

sim_uo_bstrp_out <- 
 coeftest(sim_uo_lm , sim_uo_bstrp) %>% 
 tidy(.,
      conf.int = T) %>% 
 mutate(.,
        Outcome = "Similarity  \nPerceptions",
        Group = "Ultra-Orthodox") %>% 
 filter(.,
        term == "treatment")




####Blind
sim_blind_lm <- 
 full_data %>% 
 mutate(.,
        similarity_blind_post = standardize(similarity_blind_post)) %>% 
 lm(similarity_blind_post ~ treatment + similarity_blind_pre + 
     na_similar_blind +
     grade.x + boy.x,
    data = .) 

sim_blind_bstrp <- 
 cluster.boot(sim_blind_lm, full_data$class.x, 
              parallel = T, R = 1000, 
              boot_type = "wild")

sim_blind_bstrp_out <- 
 coeftest(sim_blind_lm , sim_blind_bstrp) %>% 
 tidy(.,
      conf.int = T) %>% 
 mutate(.,
        Outcome = "Similarity  \nPerceptions",
        Group = "Visually Impaired") %>% 
 filter(.,
        term == "treatment")


#### Plot Dis-aggregated Models-----

# Combine all models together
disaggregate_models <- 
 rbind(therm_arab_bstrp_out, therm_imgrnt_bstrp_out, 
       therm_uo_bstrp_out, therm_blind_bstrp_out,
       sim_arab_bstrp_out, sim_imgrnt_bstrp_out,
       sim_uo_bstrp_out, sim_blind_bstrp_out) %>% 
 mutate(.,
        se = paste0("(", round(std.error,2),")"),
        pese = paste(round(estimate,2), se))

# set order
disaggregate_models$Group <- 
 factor(disaggregate_models$Group, 
        c("Visually Impaired",
          "Ultra-Orthodox",
          "Immigrants",
          "Arabs"))

disaggregate_models$Outcome <- 
 factor(disaggregate_models$Outcome, 
        c("Thermometers", "Similarity  \nPerceptions"))

ggplot(disaggregate_models, aes(x=estimate, y=Group,
                             label = estimate, color = "dodgerblue4")) + 
 geom_errorbar(aes(xmin=conf.low, xmax=conf.high), width=.1) +
 geom_vline(xintercept = 0,
            color = "black",
            linetype = "dashed",
            size = 0.4) +
 geom_point()+
 facet_grid(cols = vars(Outcome))+
 geom_text(aes(label= pese), 
           position = position_dodge(width = .4),
           hjust = -1.5, 
           vjust = -2.2,
           family = "Times",
           size = 2) +
 labs(x = "",
      y = "")+
 xlim(-.2,1)+
 scale_color_manual(values = c("dodgerblue4"))+
 theme(text = element_text(size = 12, family = "Times"),
       legend.position = "none",
       panel.grid.major = element_blank(),
       panel.background = element_blank(),
       panel.spacing = unit(3.5, "lines"))



##Figure A13: Interaction diagnostic-----

# Diagnose common support and linearity
interflex(estimator = "raw",
          Y = "therm_index_post", 
          D = "treatment", 
          X = "therm_index_pre", data = full_data, 
          Ylabel = "Post-Treatment Thermometer", Dlabel = "Treatment", 
          Xlabel="Pre-Treatment Thermometer", main = "",
          na.rm = T,
          theme.bw = T, show.grid = F,
          cex.main = 1.2, ncols=2)

##Figure A14: Moderating Effect of Pre-treatment Thermometor-----
# Plot interaction 
therm_int <-
interflex(Y = "therm_index_post",
                       D = "treatment", X = "therm_index_pre", 
                       Z = c("grade.x", "boy.x", "na_therm_arab",
                             "na_therm_imgrnt", "na_therm_uo", "na_therm_blind"),
                       data = full_data, 
                       estimator = "binning",
                       nboots = 1000, nsimu = 1000,
                       nbins = 3,
                       cl = "class.x",
                       vartype = "bootstrap",
                       # vcov.type = "robust",
                       na.rm = T,
                       main = "Marginal Effects")

plot(therm_int,
     theme.bw = T, show.grid = F,
     Ylabel = "Post Thermometer", 
     Dlabel = "Treatment", Xlabel = "Pre-Treatment Thermometer",
     cex.lab = .5)


##Figure A15: Pre post Arab specific outcomes-----
# Create an arab attitude index
full_data <-
 full_data %>% 
 mutate(.,
        arab_therm_pre_stn = standardize(arab_therm_pre),
        arab_therm_post_stn = standardize(arab_therm_post),
        similarity_arabs_pre_stn = standardize(similarity_arabs_pre),
        similarity_arabs_post_stn = standardize(similarity_arabs_post),
        contact_arab_pre_stn = standardize(contact_arab_pre),
        contact_arab_post_stn = standardize(contact_arab_post),
        arab_index_pre = standardize((arab_therm_pre_stn + 
                                       similarity_arabs_pre_stn+ 
                                       contact_arab_pre_stn)/3),
        arab_index_post = standardize((arab_therm_post_stn + 
                                        similarity_arabs_post_stn+ 
                                        contact_arab_post_stn)/3))

# Bind vector of pre and post arab feelings
pre_arab_index<- cbind(full_data$arab_index_pre, 
                       "pre", full_data$treatment)
post_arab_index <- cbind(full_data$arab_index_post, 
                         "post", full_data$treatment)

pre_post_arab_index <- as.data.frame(rbind(pre_arab_index, 
                                           post_arab_index)) %>% 
 rename(.,
        arab_index = V1,
        time = V2,
        treat = V3) %>% 
 mutate(.,
        arab_index = as.numeric(arab_index),
        Treatmeant = ifelse(treat ==1,"Treatment", "Control"),
        Time = case_when(
         time == "post" ~ "Post Treatment", 
         time == "pre" ~ "Pre Treatment"))


pre_post_arab_index %>% 
 filter(.,
        !is.na(arab_index)) %>% 
 group_by(.,Time, Treatmeant) %>%
 mutate(.,
        grp_mean = mean(arab_index),
        Time = case_when(
         Time == "Pre Treatment" ~ 1,
         Time == "Post Treatment" ~ 2
        )) %>% 
 ggplot(., aes(x = Time, color = Treatmeant)) +
 geom_line(aes(y=grp_mean, color=Treatmeant)) +
 geom_point(aes(y=grp_mean, color = Treatmeant, shape = Treatmeant), size = 4) +
 geom_point(aes(y = arab_index, color=Treatmeant), alpha = 0.15) +
 geom_jitter(aes(y = arab_index),  alpha = 0.15) +
 labs(x = "",
      y = "Attitudes toward Arabs Index",
      color = "",
      shape = "") +
 scale_x_continuous(breaks = 1:2,
                    limits = c(0.5,2.5),
                    labels = paste0(c("Pre", "Post"))) +
 guides(colour = guide_legend(reverse=T),
        shape = guide_legend(reverse=T)) +
 scale_color_manual(values = c("dodgerblue3", "firebrick3")) +
 theme(text = element_text(size = 16, family = "Times"),
       legend.position = "bottom",
       legend.key=element_blank(),
       panel.grid.major = element_blank(),
       panel.background = element_blank())
